4/11/24
##################
# LOAD LIBRARIES #
##################
library(tidyverse)
library(Seurat)
library(SeuratObject)
library(tidyseurat)
library(cowplot)
library(ggrepel)
library(png)
library(grid)
library(patchwork)
library(openxlsx)
source("../../bin/spatial_visualization.R")
source("../../bin/plotting_functions.R")
#########
# PATHS #
#########
input_dir <- "../../results/08_spatial_dist/"
result_dir <- "../../results/09_figures/"
if( isFALSE(dir.exists(result_dir)) ) { dir.create(result_dir,recursive = TRUE) }
epi_clus <- "^5$|^6$|^7|^9" # non-filt
#############
# LOAD DATA #
#############
# DEGs_table <- read_csv(paste0(input_dir,"subset_100/DGEs_condition_wilcox_epi_100.csv")) %>% filter(groups == "DMPA")
DEGs_table <- read_csv(paste0("../../results/05_DGE_clusters_st_data/","DGEs_clusters_wilcox.csv"))
DATA <- readRDS(paste0(input_dir,"seuratObj_spatial_dist.RDS"))
# fig.asp=5/10,
###############
# FILTER DATA #
###############
keep <- "Basal|Upper IM|Lower IM|Superficial"
names <- c("Superficial","Upper IM","Lower IM","Basal")
DAT <- DATA %>%
arrange(spatial_dist) %>%
filter( sp_annot == "epi" ) %>%
filter(!(grepl("\\d", .$layers))) %>%
mutate(groups = if_else(groups =="ctrl", "Control", .$groups))
################
# COLOR PALLET #
################
# centroids2d <- as.matrix(t(t(DAT@reductions$umap_harmony@cell.embeddings) %*% mm)/Matrix::colSums(mm))
lable_df <- DAT %>%
group_by(Clusters) %>%
select(Clusters, contains("umapharmony")) %>%
summarize_all(mean)
col1 <- set_names(c("#E41A1C","#FF7F00","#C77CFF","#984EA3"), c("6","9","7","5"))
col2 <- c("grey90","grey70", "orange3", "firebrick", "purple4")
############
# PLOTTING #
############
p1 <- plot_clusters.fun(DAT, cluster="Clusters", txt_size = 10, dot_size = 0.05,
color = col1, red = "umapharmony", lable_size = 3) + theme_void() +
xlab("UMAP 1") + ylab("UMAP 2") +
#coord_equal() +
theme(legend.position = "none",
panel.border = element_rect(colour = "black", fill=NA),
plot.title = element_blank(),
text = element_text(size = 7),
plot.margin = unit(c(5,2,1,7), "pt"), #t,r,b,l c(.1,.1,-.4,-.1)
axis.title.x.bottom = element_blank(),
axis.title.y.left = element_text(margin = margin(r = 3), angle = 90)
)
p2 <- plot_genes.fun(DAT, gene="spatial_dist", point_size = .05,
col = col2, red = "umapharmony", lable_size = 3) +
xlab("UMAP 1") + ylab("UMAP 2") +
#coord_equal() +
theme(panel.border = element_rect(colour = "black", fill=NA),
plot.title = element_blank(),
text = element_text(size = 7),
plot.margin = unit(c(1,2,1,7), "pt"), #t,r,b,l c(.1,.1,-.4,-.1)
axis.title.x.bottom = element_text(margin = margin(t = 3)),
axis.title.y.left = element_text(margin = margin(r = 3), angle = 90)
)
# dev.new(width=1.3, height=1.9, noRStudioGD = TRUE)
(A <- plot_grid(ncol = 1, #labels = c('A'),
plotlist = list(p1, p2)) )
####################
# SPATIAL DISTANCE #
####################
# dev.new(width=2.3, height=1.9, noRStudioGD = TRUE)
col1 <- set_names(c("#E41A1C","#FF7F00","#C77CFF","#984EA3"), c("6","9","7","5"))
col2 <- c("grey70", "orange3", "firebrick", "purple4")
scale_color_pseudotime <- function(palette = col2, direction = 1, ...) {
# get normalized values of spatial dist
pseudotime <- as.data.frame(list("sp"=DAT$spatial_dist))
pseudotime <- arrange(pseudotime, sp)
x <- rowMeans(pseudotime)
x <- x/max(na.omit(x))
# generate colors
pal <- palette
my_pal <- colorRampPalette(pal)(99)[x * 98 + 1]
# map them as ggplot pallet
scale_color_gradientn(colors = my_pal, na.value = "grey90", ...)
}
feat <- list("Clusters", "spatial_dist")
p <- map(feat,
~plot_spatial.fun(DATA,
geneid = .x,
sampleid = c("P118", "P097"),
zoom = "zoom",
lab = F,
save_space = F,
col = col1,
alpha = 1,
ncol = 1,
#annot_col = "#dbd9d9",
annot_line = .1,
img_alpha = 0,
point_size = .35) + #.3
#annotate("text", x = 500, y = 450, size = 4, label = "Control") + # col = "#CB251C",
theme(legend.position = "none",
panel.spacing.y = unit(-3, "lines"),
plot.margin = unit(c(-.3, -.4, -.1, -.9), "lines")) )
p[[2]] <- p[[2]] + scale_color_pseudotime(palette = col2)
# dev.new(width=2.6, height=1.9, noRStudioGD = TRUE)
B <- plot_grid(ncol = 2, #labels = c('B'), hjust = -.2,
plotlist = p)
(B <- B + draw_label("DMPA", x = c(.05), y = c(.48), hjust = 0, size = 9) +
draw_label("Control", x = c(.1), y = c( .94), hjust = 0, size = 9) +
draw_label("Gene exp. similarity", x = c(.53), y = c( .94), hjust = 0, size = 8) )
# dev.new(width=3.35, height=1.9, noRStudioGD = TRUE)
(A_B <- plot_grid( A, B, ncol=2, rel_widths = c(.55,1), labels = c('A','B'), hjust = -.2))
ggsave("../Figures/03/Fig_03AB.png", A_B, width = 3.35, height = 1.9, bg = "white", dpi = 300)
####################
# DOTPLOT FUNCTION #
####################
ord <- c("Superficial","Upper IM","Lower IM","Basal", "8","3","4","0","2","1")
clus_1 <- c("#F40014","#FF8A1E","#C77CFF","#A750B2") %>% set_names(., ord[1:4])
ID <- c("P118", "P105", "P080", "P031", "P097", "P108", "P114", "P107")
col_gr <- set_names(c("#813340", "#AC5A64", "#DA84A4", "#FFB1BA","#125E6B","#3A7D8B", "#80C2D0", "#C0DFEE"),ID)
layer_dotplot.fun <- function(DATA, feat,
facet = TRUE,
txt_size = 5,
density = FALSE,
x_max=NULL,
col){
rects <- DATA %>%
group_by(layers) %>%
summarise(., ystart=min(spatial_dist, na.rm=T), yend=max(spatial_dist, na.rm=T),
Q1=quantile(spatial_dist, probs = 0.009, na.rm=T),
Q3=quantile(spatial_dist, probs = 0.90, na.rm=T)) %>%
filter(!(is.infinite(.$ystart))) %>%
mutate(Q1 = ifelse(.$Q1 == min(.$Q1), 0,.$Q1)) %>%
mutate(Q3 = ifelse(.$Q3 == max(.$Q3), max(.$yend),.$Q3)) %>%
mutate(Q3 = ifelse(.$layers == "Lower IM", .$Q3+.3,.$Q3)) %>%
#mutate(Q1 = ifelse(.$layers == "Basal", .$Q1+.3,.$Q1)) %>%
arrange(layers) %>% ungroup()
# rects <- DATA %>%
# group_by(layers) %>%
# summarise(., ystart=min(spatial_dist, na.rm=T), yend=max(spatial_dist, na.rm=T),
# Q1=quantile(spatial_dist, probs = 0.1, na.rm=T),
# Q3=quantile(spatial_dist, probs = .9, na.rm=T)) %>%
# filter(!(is.infinite(.$ystart))) %>%
# mutate(Q1 = ifelse(.$Q1 == min(.$Q1), 0,.$Q1)) %>%
# mutate(Q3 = ifelse(.$Q3 == max(.$Q3), max(.$yend),.$Q3)) %>%
# #mutate(Q3 = ifelse(.$layers == "Basal_2", .$Q3+.3,.$Q3)) %>%
# arrange(layers) %>% ungroup() %>%
# filter(., grepl(keep, .$layers))
DAT <- DATA %>%
#filter(., grepl("SubMuc", DATA$sp_annot) & !(grepl(epi_clus, DATA$Clusters))) %>%
#filter(., (grepl(keep, DATA$layers))) %>%
mutate(., FetchData(., vars = c(feat), slot = "data") ) %>%
select(groups, ID="orig.ident", layers, all_of(c(feat)), spatial_dist) %>%
mutate(ID = factor(.$ID, levels = names(col_gr)))
dens <- ggplot() +
geom_density(data = DAT %>% filter(., .[[feat]] != 0 ), aes(y=spatial_dist, color=ID)) +
scale_y_continuous(expand = c(0, 0)) + ggpubr::clean_theme()+
scale_color_manual(values = col_gr) +
theme(panel.border = element_blank(),
plot.margin = unit(c(0,0,0,0), "inches"))
if(facet == TRUE){facets <- facet_wrap(~groups, ncol = 2)
w <- 3}
else{facets <- NULL
w <- 1}
dot <- ggplot() +
#ggtitle(feature) +
geom_rect(data = rects, alpha = 0.1, show.legend=FALSE,
aes(xmin = -Inf, xmax = Inf, ymin = Q1, ymax = Q3, fill = layers)) +
geom_jitter(data = DAT, aes(x=.data[[feat]], y=spatial_dist, col=layers),
width = 0.1, alpha = 0.7, size=.3) + #
scale_fill_manual(values = col) +
scale_colour_manual(values = col) +
coord_cartesian(clip = "off") +
guides(col = guide_legend(override.aes = list(size=2), keyheight = .7, keywidth = .7, title = "Cluster")) +
scale_y_reverse(expand = c(0, 0)) +
#scale_y_continuous(expand = c(0, 0)) +
#scale_x_continuous(expand = c(0, 0)) +
{if(!(is.null(x_max))){xlim(-.5, x_max)}} +
facets +
my_theme + ylab("Similarity in gene expression") +
theme(legend.position="top", legend.justification = "left",
#legend.position = c(.1, 1.15), legend.direction = "horizontal",
plot.background = element_rect(fill = "transparent", colour = NA),
plot.margin = unit(c(2,.0,-2,.5), "lines"),
panel.border = element_blank(),
text = element_text(size = txt_size),
axis.line = element_line(linewidth = .3),
#axis.title.y = element_blank(),
axis.text.y = element_blank(),
axis.title.x.bottom = element_text(margin = margin(t = .5)),
axis.title.y.left = element_text(margin = margin(r = 5)),
panel.grid.major = element_line(linewidth = 0.2),
panel.grid.minor = element_line(linewidth = 0.1))
if(density == TRUE){dot <- dot + dens + plot_layout(ncol=2, nrow=1, widths=c(w, 1), guides = "collect")}
return(dot)
}
###############
# FILTER DATA #
###############
keep <- "Basal|Upper IM|Lower IM|Superficial"
DAT <- DATA %>%
filter(!(grepl("\\d", .$layers))) %>%
mutate(groups = if_else(groups =="ctrl", "Control", .$groups))
##############################
# CLUS DOTPLOT PER EPI LAYER #
##############################
# dev.new(height=1.9, width=6.7, noRStudioGD = TRUE)
# dev.new(height=3.2, width=3.54, noRStudioGD = TRUE) #with legend
feat <- c("KRT78","TGM1","KRT6C","MIR205HG" )
dot_fig <- map(feat,
~layer_dotplot.fun(DAT, .x,
facet = F,
txt_size = 9,
col = clus_1)) # , x_max = 6
# dot_fig[[1]]
C_1 <- wrap_plots(dot_fig, ncol = 4, guides = "collect" ) +
plot_layout(axis_titles = "collect") &
theme(#legend.position = c(.1, 1.15), legend.direction = "horizontal",
plot.margin = unit(c(.2,0,0,.2), "lines"),
legend.position = 'top',
legend.justification = "right",
legend.background = element_rect(colour ="gray", linewidth=.2),
legend.margin=margin(1,2,1,1), # moves the legend box
legend.box.margin=margin(1,1,-4,0), # moves the legend
#legend.box.margin=margin(-30,0,-15,0),
# axis.title.y.left = element_text(margin = margin(r = 5))
)
C_1 <- plot_grid(C_1, NULL, rel_widths = c(1,.01))
ggsave("../Figures/03/Fig_03C1.png", C_1, width = 6.7, height = 2.4) #, dpi = 300
############################
# CLUS EXPRESION ON TISSUE #
############################
DAT <- DATA %>%
filter(grepl("P118", .$orig.ident))
col_feat <- c("#EFEDF5", "#DADAEB", "#BCBDDC", "#9E9AC8", "#807DBA", "#6A51A3", "#54278F", "#3F007D") # Purples
p <- map(feat,
~plot_st_feat.fun( DAT,
geneid = .x,
txt = F,
save_space = F,
zoom = "zoom",
col = col_feat,
alpha = .9,
ncol = 2,
#annot_col = "#dbd9d9",
annot_line = .1,
img_alpha = 0,
point_size = .6))
# dev.new(height=1, width=6.7/4, noRStudioGD = TRUE)
# p_list[[1]]
# dev.new(width=6.7, height=1, noRStudioGD = TRUE)
p_list <- map(p, ~ .x + theme(plot.margin = unit(c(-35, -2, -30, -15), "pt"), #t,r,b,l
#panel.background = element_rect(fill = "white"),
panel.spacing.y = unit(-1.5, "lines"),
panel.spacing.x = unit(-2, "lines"),
legend.key.height = unit(0.2, "cm"),
legend.box.margin=margin(-20,0,0,0),
legend.title = element_blank(),
legend.position = "none"
) )
#########################
# COMBINE PANEL A AND C #
#########################
# patchwork does not respect margins when nesting plots, use cowplot instead!
# dev.new(width=3.54, height=3.2, noRStudioGD = TRUE)
c <- plot_grid( plotlist = p_list, ncol = 4, rel_heights = c(1), rel_widths = c(1))
C_2 <- plot_grid(NULL, c, rel_widths = c(.05,1))
(C <- plot_grid( C_1, C_2, ncol=1, rel_heights = c(3,.9), labels = c('C'), hjust = -.2))
ggsave("../Figures/03/Fig_03C.png", C, width = 3.4, height = 3.2, bg = "white", dpi = 300)
# dev.new(width=3.4, height=1.9+3, noRStudioGD = TRUE)
A_B_C <- plot_grid( A_B, C, ncol=1, rel_heights = c(.6,1))
plot_grid( A_B_C, ncol=2, rel_widths = c(1,.03))
################################
# DOTPLOT MARKER GENE FUNCTION #
################################
library(ggnewscale)
get_order <- function(df, markers, markers_id, value){
m <- set_names(markers, markers_id)
value <- enquo(value)
df <- df %>%
mutate(marker_id = map_chr(.$marker, ~names(m)[which(m == as.character(.x))][1])) %>%
nest(., .by=c(groups)) %>%
mutate(., data = pmap(., ~filter(..2, marker_id == .[[1]]))) %>%
mutate(., data = pmap(., ~arrange(..2, desc(!!(value))))) %>%
unnest(cols = c("data")) %>%
#mutate(marker = fct_reorder2(marker, Avg, Pct))
mutate(row_id = cur_group_rows(), .by=c(groups))
return(df[["marker"]]) # use rev() to reorder the dots depending on horizontal/vertical
}
get_df_for_dotplot.fun <- function(obj, markers, groups="Clusters", gr_col, rect=TRUE, alpha=.2, gr_lvl=NULL){
col <- c("#EFEDF5", "#DADAEB", "#BCBDDC", "#9E9AC8", "#807DBA", "#6A51A3", "#54278F", "#3F007D") #
markers_lvl <- set_names(seq_along(markers), markers)
gr <- unique(obj[[groups]])[[1]]
gr_lvl <- set_names(seq_along(gr), if(is.null(gr_lvl)){gr}else{gr_lvl})
df <- obj %>%
mutate(., FetchData(., vars = c(markers)) ) %>%
as_tibble() %>%
select(., .cell, any_of(c(groups, markers))) %>%
pivot_longer(., cols = -c(".cell", groups),
names_to = "marker", values_to = "values") %>%
#mutate(Clusters = factor(.$Clusters, levels=names(clus_lvl))) %>%
group_by(!!sym(groups), marker) %>%
summarise(Avg = mean(values),
Avg_ = mean(values)/max(values), # gives the median as a percentage of max value among spots if .5 you have a normal dist.
Pct = sum(values > 0) / length(values) * 100, .groups="drop") %>%
#mutate(Pct = .$Pct*) %>% # scaling
#group_by(marker) %>%
#nest()
mutate(Avg_p = Avg/max(Avg), .by = "marker") %>%
mutate(marker = factor(.$marker, levels=get_order(., markers, markers_id, Avg))) %>%
mutate(., ymin = gr_lvl[as.character(.[[groups]])]-.5,
ymax = gr_lvl[as.character(.[[groups]])]+.5) %>%
mutate(., xmin = markers_lvl[as.character(.$marker)]-0.5,
xmax = markers_lvl[as.character(.$marker)]+0.5) %>%
mutate(!!sym(groups) := factor(unlist(.[,1]), levels=names(gr_lvl)))
# df <<- df
df_rect <- df %>% select(1, ymin, ymax) %>% slice_max(., n=1, order_by=ymin, by=groups, with_ties=F)
p <- ggplot(df, aes(x=marker, y=.data[[groups]])) +
{if(rect == TRUE)
list(new_scale_fill(),
annotate(geom = "rect", ymin=df_rect$ymin, ymax=df_rect$ymax,
xmin=rep(0, length(gr)), xmax=rep(length(markers_lvl)+1, length(gr)),
fill = gr_col[1:length(gr)], colour = NA, alpha = alpha),
coord_cartesian(expand = F) )} +
new_scale_fill() +
geom_point(aes(size = Pct, fill = Avg_, col = Avg_), shape = 21) + #"#6A2595"
scale_fill_gradientn(colours = col, aesthetics = c("fill","colour"),
guide = guide_colorbar(ticks.colour = "black",
frame.colour = "black"),
name = "Average\nexpression") +
#scale_size(max_size = 3) +
scale_size_area(max_size = 3) +
theme_bw() + labs(size="% expressed") +
guides(size = guide_legend(override.aes = list(color = "black") )) +
theme(text = element_text(size=10),
axis.text.x = element_text( angle=45, hjust=1, color="black"),
axis.text.y = element_text( color="black",hjust=.9 ),
legend.text = element_text(size=4),
legend.margin=margin(0,0,0,0), #t,r,b,l
legend.box.margin=margin(0,-6,-10,-8),
axis.title = element_blank(),
panel.grid.major = element_line(linewidth = .1),
#panel.grid.minor = element_blank(),
axis.ticks.y = element_blank(),
plot.margin = unit(c(.2,0,.2,0),units = "cm") #trbl
)
return(p)
}
###################################
# DOTPLOT TOP 10 MARKER GENES EPI #
###################################
library(ggnewscale)
# dev.new(width=6.5, height=3, noRStudioGD = TRUE)
gr_col <-c("#EB5253","#FA9F50","#D7A1FF","#984EA3")
markers <- c("RNASE7", "HMOX1", "DUOXA2", "DHRS9", "CHAC1", "PLA2G4D", "CEACAM7", "VSIG10L", "FAM3D", "LINC02487", "ENDOU", "LTB4R", "TMEM184A", "FAM83G", "ALOX15B", "COL7A1", "STMN1", "THSD4", "FGFR2", "LAMA5")
markers_id <- c(rep("Superficial",5), rep("Upper IM",5),rep("Lower IM",5),rep("Basal",5))
groups="layers"
gr_lvl <- rev(c("Superficial", "Upper IM", "Lower IM", "Basal"))
D <- DATA %>%
filter(grepl(epi_clus, .$Clusters) & sp_annot == "epi") %>%
get_df_for_dotplot.fun(., markers, groups="layers", gr_col, rect=T, gr_lvl=gr_lvl) +
theme(legend.position="bottom") + guides(fill = guide_colourbar(barwidth = 5, barheight = .5),
size = guide_legend(keyheight = 1, keywidth = .1 )
) +# horizontal
theme(#legend.box.margin=margin(-10,0,0,0),
panel.border = element_rect(colour = "gray", fill=NA),
legend.spacing.x = unit(1.5, 'pt'),
legend.text = element_text(size = 8),
legend.title = element_text(size = 8),
#legend.justification = "right",
legend.margin=margin(0,0,0,0), # moves the legend box
legend.box.margin=margin(-13,1,-5,-30), # moves the legend
axis.text.y = element_text(margin = margin(r = -1)) )
# dev.new(width=8, height=3, noRStudioGD = TRUE) # horizontal
# dev.new(width=4, height=5.5, noRStudioGD = TRUE) # vertical
# ggsave(paste0("../Figures/05/","Marker_genes_epi_h", ".pdf"), plot, width = 6, height = 3)
# ggsave(paste0("../Figures/05/","Marker_genes_epi_v", ".pdf"), plot, width = 4, height = 5.5)
# dev.new(width=3.4, height=2, noRStudioGD = TRUE) # horizontal
(D <- plot_grid(D, labels = c('D'), vjust = .9, hjust = -.2))
#############################
# COMBINE ALL FIGURE PANELS #
#############################
# dev.new(width=3.54, height=7, noRStudioGD = TRUE)
Figure3 <- plot_grid( A_B, C, D, ncol=1, rel_heights = c(.6,.85,.55))
ggsave("../Figures/03/Figure03.png", Figure3, width = 3.54, height = 6.7, bg = "white", dpi = 1000)
ggsave("../Figures/03/Figure03.pdf", Figure3, width = 3.54, height = 6.7, bg = "white", dpi = 1000)
Figure3